AN IMPROVED SYMMETRIC POSITIVE-DEFINITE FINITE 
ELEMENT METHOD FOR THE COMPLEX HELMHOLTZ 

EQUATION 

RUSSELL B. RICHINS 



Abstract. Most discretizations of the Helmholtz equation result in a system of 
linear equations that has an indefinite coefficient matrix. Much effort has been put 
into solving such systems of equations efficiently. In a previous work, the current 
author and D.C. Dobson proposed a numerical method for solving the complex 
Helmholtz equation based on the minimization variational principles developed by 
Milton, Seppecher, and Bouchitte. This method results in a system of equations 
with a symmetric positive definite coefficient matrix, but at the same time requires 
solving simultaneously for the solution and its gradient, resulting in a much larger 
system of equations. Herein is presented a method based on the saddle point vari- 
ational principles of Milton, Seppecher, and Bouchitte, which produces symmetric 
positive definite systems of equations, but eliminates the necessity of solving for 
the gradient of the solution. 

1. Introduction 

The Helmholtz equation 

V • LVu = Mu, 

is useful in modeling wave propagation in problems arising from many different physi- 
cal situations. Suppose we wish to solve the Helmholtz equation in a domain Q C M. d , 
and assume that L and M are complex valued functions. A common source of nu- 
merical methods for solving this equation is the variational principle 

(1) f [-LVu-Vv - Muv}dx = VuG^fi). 

Jn 

Since this is a stationary principle, the resulting system of equations is indefinite, 
and is generally more difficult to solve than a system of equations having a positive 
definite coefficient matrix. 

In their paper [7] , Milton, Seppecher, and Bouchitte developed variational principles 
that apply to the Helmholtz equation above, as well as the time-harmonic Maxwell 
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equations and the equations of linear elasticity in lossy materials. To derive these 
variational principles, we first define the dual variable 

v = iLVu. 

Then 

/ L \ ( Vu \ ( LVu \ ( -iv 

or equivalently, 
where 



M j \ u } \ Mu I \ -iV-v 

ZF, 



i 



JF 



Vu\ -iv \ _ L 

u r y \ -iv-v r \ o m 



For a complex quantity z, we will write z' = Re(2;) and z" = Im(^). Taking real and 
imaginary parts, the constitutive relation becomes 

Q> = Z'F - Z"T" and Q" = Z'T" + Z" F , 
which can be written in matrix form as 

Q" \ ( Z" Z' \ ( F 



( 2 ) \ q> j ~ y z' -Z" ) \ F' 

Solving this relation for the imaginary parts of J 7 and Q, we find that 

P) (£)=/:(;£, 

where 

r _fz" + z\z")- x z' z\z"^- x 



Z")- l Z' (z")- 1 
The matrix £ is positive definite as long as Z" is positive definite (see [7]). 

The previous approach in [8] was to use this constitutive relation and the correspond- 
ing energy functional 

F \ ( z" + z\z"y l z' z\z"Y l \(F\, 
-g )'\ {z")- x z' {z"Y l )\-Q'j 

to formulate a numerical method. When this variational principle is discretized by 
the finite element method, the result is a system of equations that has the block 
form 

A\ v4j Aj 
A A A 2 A\ 
A* An 












H 




v'i J 
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A similar system of equations must be solved to find approximations for u" and v'. 
In all, to find u' and u", one must solve two positive definite systems of equations of 
size 3N x 3iV, where N is the number of nodes in the computational grid. 

In this work, we develop a new method based on the saddle point variational prin- 
ciples in [7J that does not require that v be solved for in order to find u. In this 
formulation, we find u' and u" by solving two N x N symmetric positive definite 
systems of equations. Since these are real matrices, the amount of work in solving 
the two systems can be expected to be less than that required to solve the complex 
N x N system that results from before taking positivity into account, since com- 
plex multiplication is more than twice as expensive than real multiplication. 

First, in Section |2j we will derive the saddle point variational principles from [7J 
upon which our method is based. In Section [3j we will discuss the details of handling 
Dirichlet, Neumann, and Robin boundary conditions with these variational princi- 
ples. Section [4] contains the derivation of a standard bound on the error incurred 
when the Helmholtz equation is solved using a finite element method that discretizes 
the saddle-point variational principle. Section [5] outlines the numerical method and 
discusses the conditioning of the system. In Section [6j we provide several numerical 
examples, as well as numerical verification of the error bound from Section |4j 



2. The Saddle Point Variational Principle 

The derivation of the saddle point variational principle from [7j follows the same 
steps presented in the introduction for the minimization principle. The difference is 
that instead of continuing to the constitutive relation (|3]), we stop at equation ([2]). 
As usual, we assume that Z" is positive definite. From this constitutive relation, we 
define the functional 

Y(u'u") = [ ( T l V ( Z ", IV ™1 dx. 



jr" J \ Z' -Z" J \ 7" 

Let u',u" G H 1 ^) be the real and imaginary parts of a solution to the Helmholtz 
equation. Let s G Hq(Q) and define 



S 



Vs 

s 



Then we have 



Y(u' + s,u") = I ' f JC + S j • ( Z z , j ( + 5 j dx 
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Y(u',u") + 2 



r 

7" 



Z" Z' 
Z' -Z" 



s 





dx + Y(s,0). 



The integral in the line above can be rewritten as 



9. 



Q" 
Q' 



S 





dx 



n 



—v 

-V -v' 



Vs 

s 



dx 



/ \-v' • Vs - V • v's] dx= / -V • [v's] dx= -v ■ ns dS = 0. 
Jn Jq Jan 



Therefore, 



Y(u' + s, u") = Y(u', u") + I S- Z"S dx, 

Jn 



and the last term must be nonnegative, since Z" is assumed to be positive definite. 
A similar calculation yields 

Y(u', u" + s) = Y(u', u") -Is- Z"S dx. 

Jn 

This shows that (u',u") is at a saddle point of the functional Y. 

Suppose that (u',u") is a saddle point of the functional Y. Then the functional 
Q(s',s") = Y(u' + s',u" + s"), defined for all s',s" G H^(Q), should have a saddle 
point at s' — s" — 0. A necessary condition for this to happen is that the first 
variation of Q should vanish. Let 



S' 



Vs' 



Then we must have 



= 



r 
r 



and S" 



Z" Z' 
Z' -Z" 



Vs" 



S 1 
S" 



dx 



= [ [T' ■ Z"S' + J=' ■ Z'S" + F' ■ Z'S 1 - P' ■ Z"S"\ dx 
Jn 



+ 



n 
Vu" 
u" 



u 



L" 
M" 



Vs' 
s' 



Vu' 

v! 



V 
M' 



Vs" 



V 
M' 



Vs' 

j 



Vu" 

u" 



L" 
M" 



Vs 
ji 



dx 



J [Vu' ■ L"Vs' + u'M"s' + Vu' ■ L'Vs" + u'M's" + Vu"- 

L'Vs' + u"M's' - Vu" ■ L'Vs" - u" ■ M"s"} dx 

= [ [(-V • L"Vu' + M"u' -V ■ L'Vu" + M'u")s' 
Jn 
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+ (-V • L'Vu + M'u' + V • L'Vu" - M"v!')s"\ dx. 

If we rewrite the equation V • LVu = Mu in terms of real and imaginary parts, we 
arrive at 

V • [L'Vu 1 - L"Vu" + i(L'Vu" + L"Vu')] = M'u' - M"u" + i(M'u" + M"u'), 
or in other words, 

V • L'Vu - V ■ L'Vu" - M'u + M"u" = 

and 

V • L'Vu" + V • L'Vu - M'u" - M"u = 0. 

Notice that the left-hand sides of these equations are just the opposites of the ex- 
pressions multiplying s' and s" in the integral above. Since the result of the integral 
must be zero regardless of the choice of s' and s", the saddle point of Y must be a 
solution to the Helmholtz equation. 



3. Boundary Conditions 

The calculations done above show that a saddle point of Y satisfying u' = f and 
u" = f" on dfl is a solution of 

V • LVu = Mu inVl 
u = f' + if" ondQ ■ 

We can also solve the Neumann problem 

f V • LVu = Mu in Q 

y v ■ n = g' + ig" on <9fi 

Let s',s" G be arbitrary test functions. Then we have 

0= f [(V • v' - V • v')s' + (-V • v" + V • v")s"\ dx 
Jn 

= [ l-v' ■ Vs' - V ■ v's + v" ■ Vs" + V ■ v"s"] dx+ [ [sV -n- s"v" ■ n] dx 
Jn Jan 

' ! ^ ) dx + [s'v' ■ n — s"v" ■ n] dS 



Q' \ S 



n \ v J \ / J on 



f ( r \ ( z" z' \( s'\ . r r , , /; , 

— J I -p,, J ■ I z > _z" J I S" ) J i sv ' n ~ s v -n\ db. 

Therefore, in order to solve the Neumann problem, we solve the weak equation 

j[ ( jn ) • ( % ) ( %, ) dx = £ l-s'g' + s"g"] dS V s' , s" E H^Sl). 
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To solve the Robin problem 

V • LVu = Mu in Q 

u + av ■ n = g on dQ ' 

we begin with the weak form of the Neumann problem, which we will write as 

° = / n (?')'(z" -z-)( S s») dx + L{"" )'(-«""») dS 

We split the boundary condition into its real and imaginary parts as 

u' + a'v ' ■ n — a"v" ■ n — g' 
u" + a'v" ■ n + d'v ■ n = g" , 



which we can write as 



u> \ ( a 1 a" \ v'-n \ _ ( g> 
u" I 1 a " ~ a ' J V ~ v " - n \ 9" 



If the matrix in the equation above is called M , then 

v ';, n ) = -M-*(i)+ M -i( i 

—v ■ n J V / V 9 

so the weak form of the equation with Robin boundary conditions is 

-z»){s») dx -L{*") M ~ l( y*" 1 " s 

The inverse of M is 
M" 1 



—a —a \ 1 I a a 

-(a'f - (a") 2 V - a " a ' ) ~ W V a " ~ a ' 

so if we require that a' be negative, the matrix that results from discretizing the left- 
hand side will have the same block form as those that result from the other boundary 
conditions. 



4. Error Bound 
We will make the following assumptions on Z": 

,.s a. there is a constant 71 such that < 71 for all i, j, and x G f2 

b. there is 72 such that Z"{x) > 72/ for all x £ f2. 
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Define the space V = [iJ 1 (f2)] 2 , endowed with the norm 

\\(u'y)\\ v = (\\u% Hn) + \\uTm { n))^ 

Also, we will assume that Vni and Vn2 are finite dimensional subspaces of H 1 (Q), 
and Vn = Vni x Vn2 is the set in which we seek our numerical solution. 

Define an energy f(s') for s' E H 1 ^) as 

M = Iy( s ',u")+q(s',u") = lJ n (r>)-( Z z '> S' ) ( %» ) dx+Q( s >,u"), 

where Q(s',u") contains terms that arise from the enforcement of boundary condi- 
tions and any inhomogeneous terms. We will write this as 

f(s')= 1 -B(s',s')-F(s',u"), 

where 

B(s',s") = [ S'-Z"S" dx 
Jn 

and F(s',u") contains the rest of the terms. If vl is a minimizer of f(s'), then 

B(u',s')=F(s',u")W s' eH\Q). 

Therefore, we can write 

/(*') = l -B{s', s>) - F(s', u") = B(u', v!) - F(u>, u") + l -B{s' , s>) - F(s', u") 

= B(u', v!) - F(u', u") + \b(s', s') - B(u', s') 
= \b(u>, u') - F(u>, u") + l -B{u\ u') - B(u', s>) + l -B{s\ s') 
= \b(u\ v!) - F(u', u") + ^B(u' - s', u' - s') 

Suppose that u' N G Vni is such that 

f(u' N ) = min f(s'). 
seVjvi 

Then 

Bill — u' N ,u' — u' N )z = min B(u' — s',u' — s')? . 

s'ev N1 

The inequalities above imply that 

V^lWU^n) < VB(s>,s>) < Cy/^\\s'\\ H i (n) V s' E H\n), 
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and therefore, 

V72|K - u' n \\hi(q) < mill C^{\\u' - s'\\ H ifQy 

Here and in what follows, C will represent any constant that does not depend on u', 
u", or the grid spacing h. 

Let F 1 be the orthogonal projection from H 1 ^) onto Vati- Then H-FiHb^i^),.?/ 1 ^)) = 
1, where H 1 ^)) is the set of all bounded linear functions from H l (Q) to 

itself. We then take s' = F\u' to obtain the inequality 

(5) |K - u' N \\ H i {n) < C\\u' - Fxu'Whi^. 
If instead we use the energy 

f(s") = \Y{u\ s")+Q(u', s") = IJ( -Z" )(s") dx+ Q( u> > 

and perform calculations similar to those above, we obtain the bound 

(6) \\u" - u'^Wmin) < C\\u" - F 2 u"\\ H i (n) , 
where F2 is the orthogonal projection from onto Vn2- 
Combining inequalities ^ and (|6]), we find that 

II// / // "Nil 2 II' 'II 2 ill" "II 2 

\\[U —U N ,U — U N )\\y — \\U — U N \\ H i( Q ) + \\U —U N \\ H i(ty 

< C (\\u - Fxv!\\ 2 m{Fl) + \\u" - F 2 u"\\Hi {n) J = C\\{u - F x v!,u" - F 2 u")\\ 2 v 
and consequently, 

\\(u',u") - (u' N - u%)\\ v < C\\(u' - F lU ',u" - F 2 u")\\ v . 

We partition Q into subregions q, each of which can be viewed as a suitably shifted 
and rotated version of a reference element e, so that there exist affine changes of 
variables Fi(x) = Bx + xi such that Fi(e) = e\. In what follows, a hat over a function 
will denote the corresponding function defined over the reference element e obtained 
by a change of variables. 

We define the seminorm | • | s by 




and a is a multiindex. 
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From [3] we get the inequality 

(8) c~ l h s ~* \w\ afil < \w\ s < ch s -^\w\ s , ei , 

where c is a constant, w = w o F^ 1 , and | • | Sj£i denotes ^ with q in place of e. 

We now recall the following lemma from [1] : 

Lemma 1 (Bramble-Hilbert Lemma). For some region Q C R d and some integer 
k > —1, let there be given a bounded linear functional 

f : H k+1 {Q) ->• R, 

satisfying \f(u)\ < S\\u\\ H k+im) for all u G H k+1 {Vt) for some 5 independent of u. 
Suppose that f{u) = for all u G P k {pt). Then there exists a constant C , dependent 
only on Q such that 

\f(u)\<C5\u\ k+1 , ueH k+1 (Q). 

Let s G {0, 1} and fix w G H s (e). Define the functionals 

Li{u) — [u — Fiu, w] s and L 2 {u) = [u — F 2 u, w] s . 

Since 

\Lj(u)\ <\u- Fju\ s \w\ s < (\u\ s + IFj&Ulwls < (||w||i?i(e) + WFjuWhi^IwIs 
< 2||w|| H -i(g)|t(;| s < 2||w|| jH -fc+i ( g- j |w| s , 

and FjU = u for polynomial functions u in V^j (j = 1, 2), we see that the Bramble- 
Hilbert Lemma applies, and there exist constants such that 

\Li(u')\ < C\w\ s \u'\ k+1 and \L 2 (u")\ < C\w\ s \u"\ k+ i, 

as long as k is small enough so that all polynomials of degree less than or equal to 
k are contained in the span of the basis functions representing u' and u" . Taking 
w = u' — F\u' in the first inequality and w = u" — F 2 u" in the second yields 

\u — Fiu'\ s < C\u'\ k +i and \u" — F 2 u"\ s < C\u"\ k+ \. 

Assuming that h < 1 and using inequality (|8]), we see that 

\v! - F lU '\ s , ei < Ch^- S \u' - Ftu'ls < Chi- S \u'\ k+1 < Ch k - s+1 \u'\ k+t 



and 



\u" - F 2 u"\ s , ei < Ch*- S \u" - F 2 u"\ s < Chi~ s \u"\ k+1 < Ch k - s+l \u"\ k+h 
Consequently, the overall error satisfies 



(u',u") - (u' N - u" N )\\ 2 v < C\\(u' - F lU ',u" - F 2 u 



"\ 112 
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< CJ2 IW - Fm'\l ei + W - F lU '\l t + \u" - F 2 u"\le t + w - iviy 
I 

<CJ2 [h 2k+2 W\l +1 , ei + h* k \u'\t +1 , ei + h^\u'f k+1 , ei + h* k \u"\l +1 , ei ] 
i 

<Ch 2 H\u'\l +1 , u + \u'f k+1 ,u)- 

We have now proved 

Theorem 1. Under the assumptions on Z" , if the solution (u',u") G [H k+1 (fl)] 2 
and the finite element subspace used in the numerical method contains [Pk(Cl)] 2 , then 
there exists a constant C such that the error satisfies 

\\{u\u") - (u' N -u" N )\\ 2 v < Ch 2k (\u>\l +1 ,n + \u"\l + i,n), 
where h < 1 is the grid spacing. 

4.1. Regularity. For the error bound above to be useful, we must have k > 1, 
which means that u',u" G H 2 (Q), at the least. Classical regularity theory, such as 
found in [5] tells us that if L is positive definite, bounded, and has C 1 regularity, 
then u',u" G H 2 (Q). 



5. The Numerical Method 

We will examine the numerical solution of the model problem 

V • LVu = Mu inVt 



u = f + if" on Q 

The first step in solving the problem is to select a set of finite element basis functions. 
The numerical examples presented here will use a rectangular grid with bilinear basis 
functions of the form 

l_\E^A)( 1 _\y=M) if \ X - Xi \<h, \y-y k \<h 



h 

otherwise 
where h is the grid spacing. 

Regardless of how the basis is chosen, we will assume that the basis functions are 
labeled as {ifjk}k=i an d we assume that the solution has the form 

u" ) - \ % + E at&k 
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where and ip'o are auxiliary functions satisfying the boundary conditions ip' = f 
and ipQ = f" on dQ. The weak form of the Euler-Lagrange equation for the saddle 
point variational principle is 

o = J ' (% ~z" ) ( s" ) dx v S '' s " E H o(fy, 

where, as usual, 

s ' = ( v / ) and s " = ( Z" ) • 

We make the substitution above for u' and u" and let s' and s" take on the value of 
each of the basis functions. In doing so, we arrive at a system of equations which 
has the block form 

( A x Al \( a'\_fh\ 
\ A 2 -A, ){a" J \ b 2 J ' 

where A\ is positive definite. Let 

Then the entries of the blocks in the coefficient matrix satisfy 

[Al]k3 = L ( ^ ) ' ( z' -z")\ V o ) dx 

and 

[Mk3 = L ( V u ) ' ( z' -z")\ v s ) dx - 

The elements of the vector b = {pi, &2) T satisfy 

[hl]j = ~ L(n) ' \z' -z"){ V o) dx 

and 

This system of equations is of saddle point type, and therefore there is a wide array of 
numerical methods that apply (see [6]). Among the simplest is the following, based 
on a Schur complement. By using this approach, we reduce the problem from solving 
an indefinite 2N x 2N system to solving two N x N positive definite systems. First, 
we rearrange the second equation, finding 

Aid' = -b 2 + A 2 a', 
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then multiply by A 1 and substitute the result into the first equation to get 

(A 1 + AlA^A 2 )a' = 61 + AlA^b 2 . 

Because A x is positive definite, the coefficient matrices in both these systems of 
equations are positive definite. 

While the matrix A 1 + A\A\ X A 2 is positive definite, it is also costly to store and to 
compute. For this reason, we use the preconditioned conjugate gradient method to 
compute the solution to the system with this coefficient matrix, since this method 
only requires the ability to perform matrix-vector multiplication with the coefficient 
matrix. As a preconditioner for A\ + A\A~{ X A 2 , we use the matrix A±. 

Systems with coefficient matrix A x appear explicitly in the algorithm, but must also 
be solved at each step when PCG is applied to the matrix A\ + A\A\ l A 2 . The 
systems of equations of the form A\x = b can also be solved using PCG. An effective 
strategy for preconditioning A\ is to write it as a sum two matrices, one arising from 
interactions of gradients of basis functions, and the other arising from interactions of 
the basis functions with themselves, then use an incomplete Cholesky factorization of 
the first of the matrices in the sum to precondition. In other words, if we write 

A 1 = Pi + P 2 , 

where 

[Pi]kj = I W> fc • L"V^j dx, 
Jn 

we use an incomplete Cholesky factorization of P\ as a preconditioner for A±. 
The algorithm is as follows: 

(1) Form the matrices Ai and A 2 . 

(2) Compute the right-hand side vectors b\ and b 2 . 

(3) Compute Wi — b± + A\A~{ x b 2 . 

(4) Solve (Ai+A^A^A^a' = w 1 using PCG with the preconditioner A 1 = R T R. 

(5) Compute w 2 = —b 2 + A 2 oi '. 

(6) Solve A\a" = w 2 by PCG with an incomplete Cholesky factorization of Pi as 
preconditioner. 

In the algorithm above, whenever it is necessary to solve a system with coefficient 
matrix Ai, we use PCG with preconditioner P\. 

This algorithm is completely implicit, and therefore is well suited for large-scale prob- 
lems. Because all that is required are sparse matrix-vector multiplications, parallel 
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implementations of this algorithm should produce a significant speedup. In partic- 
ular, this algorithm could be implemented on a GPU cluster, where many graphics 
processing units (each of which contains many processors) are used in parallel to 
perform very fast computations. 

If a completely implicit algorithm is not necessary, we may perform a Cholesky 
factorization on A\. Then systems with A\ as the coefficient matrix can be solved 
easily by performing back substitution twice. 



5.1. Conditioning. The primary goal of this section is to determine how the con- 
ditioning of the matrix A depends on the coefficients L and M in the Helmholtz 
equation. To better understand this relationship, we wish to see how the choices of 
L and M affect the eigenvalues of the matrix 

Z" Z' 
Z' -Z" 



(9) 



As we have seen, we may factor a 2 x 2 block matrix as 

A — BD~ X C 





B > 










i- a 





D~ l C 



to conclude that 



det ( C D ) = det ( D ) det ( A " BD- l C). 



Assuming det(Z') ^ 0, this implies that 



z" - \i z' \ , z' -z"-\i 



det ^ z , _ z „ _ X] . J - ( 1) det ^ z „ _ M z , 
= (-l) d+1 det(Z / )det(Z / - (-Z" - XI)(Z')-\Z" - XI)) 

= (-l) d+1 det(Z / )det(Z / + z"(z')- l z" - XZ" {z'y 1 + x{z'y l z" - x 2 (z')- 1 )- 

If we write Z = diag(ci, . . . , Cd+i), then the eigenvalues of ^ must satisfy 

-AVy)- 1 + 4 + WOVi)" 1 = o, 

and therefore, 

A = ±^-) 2 + (^) 2 = ±N. 

This means that in order to have a well conditioned system, we want to ensure that 
the entries of Z are large in terms of absolute value, or equivalently, we want \L\ and 
\M\ to be large. 



Regardless of how L and M are chosen, we are still discretizing a differential operator, 
and as such, the resulting matrix can be expected to suffer from poor conditioning. 
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Figure 1. The distribution of the eigenvalues of A\ + A\A~{ X A 2 (left) 
and the real parts of the eigenvalues of A[ X {A\ + A^A^ 1 A 2 ) (right) for 
an example with 30 2 nodes in the computational grid. 

In the numerical algorithm outlined above, we suggest that A\ be used as a precon- 
ditioner. In Figure hi we see the distribution of the eigenvalues of A\ + A\A\ X A 2 
and A^ l (A 1 + A\A~{ A 2 ) for an example where the real and imaginary parts of L 
and M take on random values in the range (0, 10). 

When we employ this method of preconditioning, the number of iterations required 
to solve a problem to a given tolerance on the size of the relative residual remains 
the same, regardless of the number of nodes in the computational grid (see Fig- 
ure gj. 

So far, we have assumed that Z" is positive definite, but it is often possible to use 
this method even when L and M do not have positive imaginary parts. A solution 
of the equation 

V • LVu = Mu 

is also a solution to 

V • e ie LVu = e id Mu, 

where 9 is a constant. Therefore, to assure that the imaginary part of Z" is positive 
definite, we can rotate them in the complex plane to assure that the new coefficients 
e L and e l9 M have positive imaginary parts. The necessary conditions on L and M 
for the method to apply are that their values lie within one half-plane. That half 
plane may then be rotated so that it becomes the upper half-plane. 

Care must be taken with solving the Neumann and Robin problems when this method 
of rotation is used, to ensure that the correct boundary conditions are enforced. For 
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Figure 2. The number of PCG iterations needed to solve a problem 
with L = 2 + .003z and M = —3 + .0004? as a function of the grid size 
N for different tolerances on the relative residual. 

example, if one desires to solve the Neumann problem 

f V • LVu = Mu in Q 
1 v ■ n = g on <9f2 ' 

the rotated version of the problem is 

f V • e ie LVu = e i9 Mu in Q 
1 v • n = e %e g on dVL ' 

where w = ie ld LWu. In Figure [3j the result of rotation on an example with L = 3 + 2i 
and M = 1 + 4z is shown. The error remains nearly constant until 9 is such that one 
of the imaginary parts of the rotated coefficients approaches zero. 



6. Numerical Examples 



In this section, we present several numerical examples solved by this method. In 
Figure |4| we solve the Helmholtz equation in a layered material, where L = 3 + 
2i, M = 1 + Ai in one material, and L = 0.5 + O.OOlz, M = 3 + 7i in the other. 
The equation is solved with Dirichlet boundary conditions provided by the auxiliary 
functions 

ip' = cos(1.5x) cos(1.5x), ipQ = sin(x) sin(y). 
In Figure [5j M has the constant value 63.9923 + 0.7039i, while L has the value 
— 0.5000 + 0.0027? in a bar of width 1/4 angled from the upper left to the lower right. 
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Figure 3. The effect of rotation on a sample problem with constant 
coefficients. The error remains nearly the same for all values of 9 that 
keep L and M in the upper half plane. 





Figure 4. The real (left) and imaginary (right) parts of a solution to 
the Helmholtz equation in a layered material. 

The equation is solved with Robin boundary conditions with a = — 1 + (l/3)i and 
g = 3.333z. 

Unfortunately, this method suffers from the usual difficulties when trying to solve 
problems with a high frequency. For example, when 

J 2 

L = -p , M = — , and u = P, 

K 
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Figure 5. The real (left) and imaginary (right) parts of a solution to 
the Helmholtz equation in a two-phase material with Robin boundary 
conditions. 



N 


h 


\\(u'-u' N ,u" 




30 


0.0345 


1.0402 x 


10" 


-4 


40 


0.0256 


5.7960 x 


10 


-5 


50 


0.0204 


3.6737 x 


10 


-5 


60 


0.0169 


2.5327 x 


10 


-5 


70 


0.0145 


1.8451 x 


10 


-5 


80 


0.0127 


1.4563 x 


10' 


-5 


90 


0.0112 


1.1154 x 


10' 


-5 


100 


0.0101 


9.0364 x 


10 


-6 



Table 1 . The error in the finite element solution for various values of 
the grid size N. 



the Helmholtz equation becomes the acoustic equation, where p is the material den- 
sity, k is the bulk modulus, u is the frequency, and P is the pressure. In figure [6j we 
solve the acoustic equation with 

p = 2 + 2i and k = 1 — 3i 

and let u range between 1 and 30, and find that the error grows rapidly with u, even 
as the ratio of uj to the grid spacing is held approximately constant. 

In Table [TJ we see the decrease in the error for an example problem as the number 
of grid points is increased. 
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omegs 

Figure 6. The error as a function of u in solving the acoustic equa- 
tion, with the ratio of u to the grid spacing held approximately con- 
stant. 

7. Conclusion 

By formulating the finite element problem through the saddle-point variational prin- 
ciples of Milton, Seppecher, and Bouchitte, we are able to solve boundary value 
problems for the complex Helmholtz equation by solving symmetric positive defi- 
nite systems of equations. The method is based on using elimination on the block 
structure of the finite element matrix to produce two smaller systems of equations, 
both of which have positive definite coefficient matrices. The systems can be solved 
using purely iterative methods, or by a combination of an iterative method with a 
Cholesky factorization. This is a marked improvement over the previous method 
based on minimization variational principles that required that the u and v variables 
be solved for simultaneously. 

This method applies to a large class of problems, especially in light of the ability to 
"rotate" the coefficients of a given problem to fit the assumptions of the algorithm. 
The symmetric positive definite nature of the matrices resulting from this method 
allow for a simple and efficient solution procedure. 

It should be emhasized that the method developed here does not only apply to the 
Helmholtz equation. In [7j, there are similar variational principles given for the time- 
harmonic Maxwell equations and the equations of linear elasticity in lossy materials. 
The ideas presented here can easily be adapted to these situations. Also, the original 



AN IMPROVED METHOD FOR THE HELMHOLTZ EQUATION 



19 



variational principles of this type, developed by Cherkaev and Gibianski in [3], can 
be used to apply this numerical method to the complex Poisson equation. 

As with the previous minimization-based method, the variational principles upon 
which this method is based remain valid as long as L and M have positive imaginary 
part, but the conditioning of the system deteriorates and the error incurred increases 
as L and M come closer to violating this condition. 

There is still more study necessary to determine the exact conditions under which 
this approach outperforms other methods already in use. Also, it is worthwhile to 
consider other boundary conditions in addition to the simple ones presented herein, 
such as a PML j9]. A forthcoming paper will address the application of this method 
to problems with a non-local boundary condition, such as those presented in [2]. 
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